Network synchronization landscape reveals 
compensatory structures, quantization, and 
the positive effect of negative interactions 

Takashi Nishikawa * ^ and Adilson E. Motter * * 

* Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA, Department of Mathematics, Clarkson University, Potsdam, NY 13699, USA, and 

* Northwestern Institute on Complex Systems, Northwestern University, Evanston, IL 60208, USA 

Published in PNAS 107, 10342-10347 (2010) 



Synchronization, in which individual dynamical units keep in pace 
with each other in a decentralized fashion, depends both on the dy- 
namical units and on the properties of the interaction network. Yet, 
the role played by the network has resisted comprehensive character- 
ization within the prevailing paradigm that interactions facilitating 
pair-wise synchronization also facilitate collective synchronization. 
Here we challenge this paradigm and show that networks with best 
complete synchronization, least coupling cost, and maximum dynam- 
ical robustness, have arbitrary complexity but quantized total interac- 
tion strength, which constrains the allowed number of connections. 
It stems from this characterization that negative interactions as well 
as link removals can be used to systematically improve and optimize 
synchronization properties in both directed and undirected networks. 
These results extend the recently discovered compensatory pertur- 
bations in metabolic networks to the realm of oscillator networks and 
demonstrate why "less can be more" in network synchronization. 

complex networks | nonlinear dynamics j stability analysis j network interac- 
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Flocking animals [l] , self-coordinated moving sensors [3] , 
bursting neurons [4], pace-matching chemical oscillators 
[5], and frequency-locked power generators [6] are some of 
many physical manifestations of spontaneous synchronization. 
Like other forms of collective phenomena 71, synchronization 
depends critically on the properties of the interaction network 
[5]. Common wisdom suggests that synchronization is gener- 
ally easier to achieve with more interactions, that synchroniza- 
tion properties change monotonically as the number of avail- 
able interactions is varied, and that certain network structures 
facilitate while others inhibit synchronization. These three ex- 
pectations, however, are all false because they ignore the pos- 
sibility of compensatory structural effects. For example, re- 
moving a link from a globally connected network makes it less 
likely to synchronize, but targeted removal of additional links 
can enhance synchronization [9l [lOl [11] [I2l [H [H [15] . Het- 
erogeneous distribution of coupling strengths or connectivity 
generally inhibits synchronization |161I17) . but when combined 
they can compensate for each other [181119] . Bringing this ar- 
gument one step further, while previous studies have focused 
mainly on positive interactions (but see Refs. [201 1211 122] ) 
— presumably because negative interactions alone generally 
do not support synchrony — it is actually easy to provide 
examples in which negative interactions help stabilize syn- 
chronous states. This is illustrated in Fig. [TJ where the net- 
work composed of black and blue links is not optimal for syn- 
chronization but the removal of the blue interactions or, alter- 
natively, the addition of interactions with negative strengths 
(red links) makes it optimal; the same is achieved by weight- 
ing the strengths of all three input interactions of each purple 
node by a factor of 2/3. However, the counter-intuitive prop- 
erties that start to emerge from such case studies currently 
lack a common in-depth explanation that is both comprehen- 
sive and predictive in nature. 




Fig. 1. Examples of compensatory network structures. Nodes represent identical 
dynamical units and links diffusive couplings. The network consisting of the black 
and blue links, all having strength is not an optimal network for synchronization, 
namely one that synchronizes over the widest range of global coupling strength (to be 
formalized in Results section). However, it can be made optimal either by 1) removing 
the blue links, 2) adding the red links with negative strength —1, or 3) scaling the 
strengths of all in-links of the purple nodes by a factor of 2/3. 

Here we show that these and other apparently disparate 
properties follow from the discovery we present below that 
networks optimal for synchronization have a quantized num- 
ber of links, in multiples of a constant that depends only on 
the number of nodes and on the connection strengths. We 
derive our results for the local stability of synchronous states 
in networks of identical units and provide evidence that our 
main results remain valid for networks of non-identical units. 
We choose to focus on optimal networks because, as we show, 
this class is very rich, can be dealt with analytically, and forms 
a multi-cusped synchronization landscape, which underlies all 
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synchronizable networks and from which suboptimal networks 
can be studied using perturbation and numerical methods. An 
immediate consequence of our quantization result is that the 
addition of links to an optimal network generally results in 
a suboptimal network, providing systematic examples of link 
removals that enhance synchronizability. Similar quantization 
results hold true for optimal networks with negative interac- 
tions, which we derive using a generalized complement trans- 
formation that maps them into networks with only positive in- 
teractions. We show that negative interactions can reveal an- 
tagonistic structural relations and counterbalance structural 
heterogeneities, with potential implications for inhibitory in- 
teractions in neural [231 1241 [25] . power-grid [6l[26|, and cell- 
regulatory networks [27] . The interactions between power-grid 
generators, for example, can be positive or negative depending 
on the inductive versus capacitive nature of the corresponding 
network elements (e.g., for the Western U.S. power grid they 
are split 97% versus 3% [2S]). 



Results 

Optimal networks for synchronization. We represent the struc- 
ture of a network with n nodes using its adjacency matrix 
A = {Aij)i<ij<:n, where Aij is the strength of the link from 
node j to node i. We consider the network dynamics governed 
by 



= F(xO+eX^A,,[H(x,)-H(x. 



[1] 



where Xi is the state vector of the ith dynamical unit, F rep- 
resents the dynamics of an isolated unit, H(xj) is the signal 
that the jth unit sends to other units [28], and e — e/d is 
the global coupling strength e > normalized by the the av- 
erage couphng strength per node d := ^ij ■ As 
a result of this normalization, the dynamics for a given e is 
invariant under scaling of ^4 by a constant, which does not 
change the network structure. This system has served as a 
workforce model to study synchronization because it allows 
analysis of the network influence without detailed specifica- 
tion of the properties of the dynamical units. For example, 
using a stability function A(-) that is independent of the net- 
work structure [T] 1281 1291 130] , the condition for a synchronous 
state xi(t) — ■■■ = x„(i) = s{t) to be linearly stable is 
A(eAi) < for i = 2, ...,n, where A2,...,An are the non- 
identically zero eigenvalues of the Laplacian matrix L defined 
by Lij = Sij J, Aik — Aij (see Materials and Methods). Thus, 
the smaller the normalized spread of the eigenvalues in the 
complex plane, which we quantify using 
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|Ai - A|^, where A := ^ Ai, [2] 



the more synchronizable the network will generally be. An- 
other measure of synchronizability is the coupling cost K : — 
£s X^i X^jT^i = Es ■ m at the synchronization threshold 
e = Es, whose minimization is equivalent to the condition 



Aa 



A3 



A„ = A > 0, 



which is also equivalent to the condition a = 0. This condi- 
tion constrains all the eigenvalues to be real, even though the 
networks can be directed and directed networks have complex 
eigenvalues in general. Condition [3] is also equivalent to the 
maximization of the range of e that allows for stable synchro- 
nization as well as to the maximization of dynamical robust- 
ness, in that it can achieve the fastest exponential convergence 




K 100 



Fig. 2. Quantization in the synchronization and cost landscapes. (A) Estimated 
by simulated annealing for networks with n = 10-12 nodes, the minimum value 
fJniin (blue dots) of a (Eq. [2j shows cusp points at periodic values of the total 
interaction strength, TTi = qj^ (Eq. [4). In this representation as a function of m, 
the results for networks with positive {Aij = 0, 1) and mixed positive/negative 
[Aij = 0, dil) interactions are identical (blue dots) and coincide with the the- 
oretical prediction (black curves) (Eq. [6j. The results for the case Aij = ibl 
(red dots) shows a similar but different quantization behavior, following Eq.[8] (B) 
Synchronization landscape Cj^iij^ as a function of both m and n, computed from 
Eq. \Q\ Through simulated annealing, a network with initially high tj (top red dot) 
evolves to a network satisfying Eq. [5]and <J = cr^iin 

(bottom red dot). (C) The 
coupling cost as a function of m and n in units of K2 (the cost for a networi< of 
two coupled units). The cost takes the lower values along the lines of discontinuity 
[m = ijfc). 



to synchronization (see Materials and Methods). (This may 
relate, for example, to the finding that heterogeneous distri- 
bution of links, which tend to make the distribution of Ai's 
heterogeneous, leads to longer transients to synchronization 
in ecological network models with strong couplings [31].) We 
emphasize that the equivalence of these conditions holds for a 
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wide variety of possible stability functions, including those for 
which the region defined by A(-) < may be finite or semi- 
infinite, or have multiple disconnected components. Having 
the maximum range, minimum coupling cost, and maximum 
dynamical robustness, the networks satisfying Eq. [3]are called 
the optimal networks for synchronization. Similar results can 
be derived for networks of non-identical units, in which the 
functions F and H are possibly different for different nodes, 
and this more general case will be discussed below. 

Quantized number of links. A surprising consequence of hav- 
ing these optimal synchronization properties is the quanti- 
zation of the number of links in the networks. We find, for 
example, that for binary interaction networks (i.e., Aij — 0, 1) 
satisfying condition [31 the number of links m = "^^i "llj^i 
is quantized to multiples of n — 1. That is, 

m = Qk :— k{n — 1), where = 1, 2, . . . , n. [4] 
This follows from the identity m — - La = X]r=2 ~ 
(n — 1)A, combined with the fact that condition [3] constrains 
the real eigenvalue A further to be an integer for networks with 
integer-valued interactions (see Supporting Information, Sec- 
tion 1). Consequently, any network with m strictly between 
these quantized values must have a > 0, and hence cannot be 
optimal. What is then the minimum a for all such networks 
with a given m? We denote this minimum by cr = aminim). 
Based on our analysis of all networks with n < 6, we conjec- 
ture that the condition to achieve this minimum is that the 
Laplacian eigenvalues (counting multiplicity) have the form 

0,fc~^^,'fc-H,.^,A:-hl, [5] 
where k is the integer satisfying < m < qk+i- Note that, 
analogously to Eq. |4]for m = qu, this condition asserts that 
the eigenvalues are not only real but also integer for any m. 
This leads to our prediction that 

o-mm(m) = ^^^ ^/ (m - qk){qk+i - m), [6] 

which is expected to be valid for binary interaction networks 
with arbitrary number of nodes and links. Indeed, Fig. [2l4 
shows that for 10 < n < 12, a simulated annealing procedure 
identified networks (blue dots) with a within 10~^ of (but 
not smaller than) (Jinin(m) predicted by Eq. [6] (black curve). 
Starting from the (optimal) fully connected network [with 
m = n(n— 1)] at the right end of the curve (Jmin(m), any initial 
link deletion necessarily makes synchronization more difficult. 
Further deletions, however, can make it easier, bringing the 
networks back to optimal periodically as a function of m at 
the cusp points and eventually leading to (optimal) directed 
spanning trees with m = n — 1 (see Supporting Video). The 
optimization of synchronization at these cusp points is analo- 
gous to the optimization of the dynamical range of excitable 
networks at the critical points of phase transitions [32]. Sim- 
ilar cusps are also observed for the Laplacian eigenvalues of 
a structurally perturbed optimal network as a function of the 
perturbation (see Supporting Information, Section 2). Note 
that, although the cost K — Ss ■ m generally depends on m, 
it is actually independent of m for optimal networks of given 
size n because the synchronization threshold Es = £s(m) com- 
pensates for any change in m. 

To reveal the intricate dependence of synchronization on 
the number of nodes and links, we now consider (Tmin as a 
function of both m and n (Fig. [2]B) based on our predic- 
tion [6] Each point on this synchronization landscape repre- 
sents all networks with a = (Tmin for a given m and n, and 
the number of such networks is expected to grow combinatori- 
ally with n (see Supporting Information, Section 3). All other 



networks with the same m and n are represented by points 
directly above that point. The evolution of one such network 
by rewiring links under pressure to optimize synchronizability 
can be visualized as a vertical descent from a point above the 
landscape toward a point on the landscape that minimizes a 
(red arrow in Fig. [2]B). The sharp valleys are observed along 
the lines m — qk, and therefore correspond to the cusp points 
in Fig. [214. Because two adjacent points on the landscape 
may include networks that do not have similar structures, it 
is surprising that one can actually move from one point to the 
next often by a simple structural modification (see Supporting 
Video). Indeed, moving in the direction of the m-coordinate 
can be achieved by simply adding or removing a link that in- 
duces the smallest increase or largest decrease in a. Along 
the line m = gfe, an optimal network can be "grown" and 
kept optimal, by adding a node and connecting any k existing 
nodes to the new node (see Supporting Information, Section 
3). The flexibility of choosing new links in this construction 
strongly suggests that optimal networks can have arbitrarily 
complex connectivity structures as the network size grows. 
Another interesting landscape appears when we compute the 
cost if as a function of m and n (Fig. [2]C) based on condi- 
tion [5] In this landscape, optimal networks lie along the lines 
of discontinuity, resulting from the discontinuous change in Sa 
that occurs when m changes from — 1 to qk and defines a 
saw-tooth function along the m-coordinate. Note that for op- 
timal networks, the cost K is independent of m, as mentioned 
above, but increases linearly with n and can be expressed as 
K = K2{n — 1), where K2 is the coupling cost for a net- 
work of two coupled units. A different but potentially related 
structure is the roughness of the time horizon considered in 
the synchronization of parallel processing units in distributed 
computing 33, 34 . 

While the presentation above focused on networks of iden- 
tical units, the main results also hold true for networks of 
non-identical units. Adopting the framework of Ref. [35] and 
developing further analysis for networks of one-dimensional 
maps, we show in Supporting Information (Section 4) that 
complete synchronization is possible even for non-identical 
units. Moreover, we show that this is possible only for net- 
works satisfying Eq. [3] in addition to the condition that the 
Laplacian matrix L is diagonalizable. Since any such networks 
must also exhibit the same quantization expressed in Eq. [4] 
we also expect cusps similar to those shown in Fig. [3 For 
each quantized number of links qk, we can show that there is 
a network that can synchronize completely despite the hetero- 
geneity of the dynamical units. Therefore, for both identical 
and non-identical units, condition [5] can be used to systemat- 
ically construct examples of suboptimal networks that can be 
made optimal by either adding or removing links. 

Stabilizing effect of negative interactions. Interestingly, the 
exact same quantization effect described above for binary in- 
teraction networks is also observed when we allow for nega- 
tive interactions and interpret m = '^jy^i ^ 
number of links. To see this, we use a generalization of the 
complement network, which we define for a given constant a 
to be the network with adjacency matrix A'^ given by 

:=(«-A,,)(l-5,,). [7] 

(This includes the special case a = 1, which for undirected 
unweighted networks corresponds to the previously studied 
case of complement graphs [36].) The transformation from 
a network to its complement maps m to an{n — 1) — m, Xi 
to an — Xi, a to — , , and thus an optimal network 
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to another optimal network when a > (see Materi- 

als and Methods). This also establishes a mapping between 
networks capable of complete synchronization of non-identical 
units, since the Laplacian matrix for such a network remains 
diagonalizable under the transformation (see Supporting In- 
formation, Section 4). The condition on a avoids (nonsyn- 
chronizable) networks having eigenvalues with negative real 
part as an artifact of the transformation. We argue that this 
generalized complement transformation is a powerful tool in 
analyzing networks with negative interactions because it re- 
duces problems involving negative interactions to those involv- 
ing only positive interactions when we choose a > max Aij. 

As an example of quantization with negative interactions, 
we consider the class of networks for which Aij = 0, ±1 and 
assume that A2 , . • . , A„ have positive real parts to ensure that 
the network can synchronize. Mapping these networks under 
the complement transformation with a = 1, we obtain pos- 
itive interaction networks with Aij — 0,1,2. Conjecture [5l 
applied to the resulting networks then leads to a prediction 
identical to Eq. [H] which we validated by simulated annealing 
for networks of size up to 12 (blue dots. Fig. [214) • Thus, 
all the results discussed above based on Eq. [S] including 
Eq. 13] and the shape of the synchronization landscape, are 
predicted to be valid even when negative interactions are al- 
lowed. Whether the networks with a — amin{m) actually do 
have negative interactions is not a priori clear because the 
synchronization landscape can be built entirely by using the 
subset of networks with only positive interactions. Our simu- 
lated annealing shows, however, that many optimal networks 
(i.e., with cr = 0) have negative interactions, as illustrated by 
an example in Fig. [T] In addition, complete synchronization 
of non-identical units is possible in the presence of negative 
interactions (see Supporting Information, Section 4). These 
networks provide clear and systematic examples of negative 
interactions that improve synchronization, as removing nega- 
tive interactions from any such network would in general push 
m off the quantized values (m = qk) and make the network 
suboptimal. 

The quantization of m and the shape of the synchroniza- 
tion landscape, though they were identical for the two exam- 
ples above, do depend critically on the constraints imposed on 
the interactions. Consider, for example, the fully connected 
networks with interaction strengths ±1. It can be shown that 
the implied constraint Aij 7^ leads to a different quantiza- 
tion compared to Eq. [6] which involves multiples of 2{n — 1) 
rather than n — 1, 



Crmm(m) = l^tj V ("^ ^ 9n-2(fe + l))('7n-2fc " m) [8] 

for g„_2(fc+i) < m < g„-2fc (red dots. Fig. [2l4). The syn- 
chronization landscape can be drastically different in some 
extreme cases. On the one hand, within the widely studied 
class of undirected unweighted networks (corresponding to a 
stronger set of constraints, Aij — Aji, Aij = 0, 1), no network 
with m < n(n— 1) satisfies the optimality condition [5] |37l [TH] . 
which indicates that aminim) > for all m < n(n — 1). On the 
other hand, within the class of weighted networks correspond- 
ing to having no constraint on the interaction strengths, an 
optimal network can be constructed for any number of links 
> n — 1 (e.g., the hierarchical networks in Ref. [10]), which 
implies that cTmin is identically zero in this particular case. 

To demonstrate that the synchronization enhancement by 
negative interactions goes much beyond the realm of opti- 
mal networks, we propose a simple algorithm for assigning 
strength —1 to directional links in an arbitrary network with 
all link strengths initially equal to +1. Our strategy is based 
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Fig. 3. Improving synchronization with negative interactions. (A) Strength — lis 
preferentially assigned to directional links (red arrows) connecting to high in-degree 
nodes (those in the shaded area) in an initially undirected network. The black links 
have strength -t-1. [B] The initial spread of the Laplacian eigenvalues is significantly 
reduced as a result (from blue to red dots). The error bars indicate the mean and one 
standard deviation of the corresponding distributions. (C) The mean and the stan- 
dard deviation (red error bars) of the percentage reduction in cr over 20 realizations 
(blue dots) of random undirected scale-free networks with 1,000 nodes, minimum 
degree 5, and the scaling exponent 7. The networks were generated by connecting 
nodes randomly according to a power-law degree distribution P{k) ~ . 

on the observation that the in-degree distribution is a main 
factor determining synchronizability [17lll8l[T9ll38ll39] . where 
the in-degree (or the total input strength) of node i is defined 
to be "^j^i Aij . Since heterogeneity in the in-degree distribu- 
tion tends to inhibit synchronization [171 1181 [T9l 1381 140) . here 
we use negative interactions to compensate for positive in- 
teractions and to homogenize the in-degree distribution. For 
a given network, we first choose randomly a node with the 
smallest in-degree and change the strength of each out-link of 
that node to —1, unless it makes the in-degree of the target 
node smaller than the mean in-degree of the original network. 
We keep strength -1-1 for the other out-links, as well as all 
the in-links. Having treated this node, we repeat this process 
for the subnetwork of untreated nodes, considering links and 
hence degrees only within that subnetwork. We continue this 
process until all nodes are treated. Applying this algorithm to 
the network in Fig. [3j4, we see that the high in-degree nodes of 
the initial network (those in the shaded area) receive all of the 
compensatory negative interactions (red arrows), reducing a 
by nearly 35% (Fig. [3]B). The effectiveness of the method was 
further validated (Fig. [3]C) using random scale-free networks 
[41] generated by the standard configuration model [42] , where 
we see more dramatic effect for more heterogeneous networks. 
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reducing a by as much as 85%. The synchronization enhance- 
ment was indeed accompanied by the homogenization of the 
in-degree distribution (see Supporting Information, Section 
5). The use of negative directional interactions in our algo- 
rithm suggests that the enhancement is partly due to link di- 
rectionality (see Ref. [13] for an enhancement method purely 
based on link directionality), which generally plays an impor- 
tant role in synchronization (see for example Refs. [21) and 
[44]). However, negative strength of interactions alone can 
also produce similar enhancement in random scale-free net- 
works when they are assigned to bidirectional links between 
hubs (see Supporting Information, Section 6). 



Conclusions 

Even though negative interactions and link removals by them- 
selves tend to destabilize synchronous states, we have shown 
that they can compensate for other instabilities, such as those 
resulting from a "forbidden" number of interactions or from a 
heterogeneous in-degree distribution. This establishes an un- 
expected relation between network synchronization and recent 
work on metabolic networks, where locally deleterious pertur- 
bations have been found to generate similar compensatory 
effects that are globally beneficial in the presence of other 
deleterious perturbations |45l I46| . For example, the removal 
of a metabolic reaction — or, equivalently, of its enzyme- 
coding gene(s) — can often be partially or totally compen- 
sated by the targeted removal of a second metabolic reaction. 
That is, the inactivation of specific metabolic reactions can 
improve the performance of defective or sub-optimally oper- 
ating metabolic networks, and in some cases it can even rescue 
cells that would be nonviable otherwise [33]. Other research 
has shown that similar reaction inactivation occurs sponta- 
neously for cells evolved to optimize specific metabolic ob- 
jectives 07] • These apparently disparate examples share the 
common property that the collective dynamics is controlled, 
and in fact enhanced, by constraimng rather than augmenting 
the underlying network. Such network-based control, includ- 
ing the conditionally positive impact of otherwise detrimental 
interactions and constraints, is most likely not unique to these 
contexts. 

In neuronal networks, for example, inhibitory interac- 
tions have been predicted to facilitate synchronous bursting 
^3, 24, 23 HH] IMl Even more compelling, it has been 

established for both animal :51 and human [521 153[ cerebral 
cortices that the density of synapses first increases and then 
decreases during early brain development, suggesting a pos- 
itive role played by link removal also in neuronal networks. 
More generally, we expect that the eigenvalue spectrum analy- 
sis that lead to our conclusions will help investigate analogous 
behavior in yet other classes of network processes governed by 
spectral properties, such as epidemic spreading and cas- 
cading failures [551 [56] . 

Taken together, the highly structured characteristics re- 
vealed by the synchronization landscape explain why the char- 
acterization of the network properties that govern synchro- 
nization has been so elusive. Numerous previous studies per- 
formed under comparable conditions have sometimes led to 
apparently confiicting conclusions about the role played by 
specific network structural properties such as randomness, be- 
tweenness centrality, and degree distribution [3^ . Our results 
indicate that at least part of these disagreements can be at- 
tributed to the sensitive dependence of the synchronization 
properties on the specific combination of nodes and links, as 
clearly illustrated by the non-monotonic, periodic structure of 
cusps exhibited by the synchronization landscape. We suggest 
that insights provided by these findings will illuminate the de- 



sign principles and evolution mechanisms of both natural and 
engineered networks in which synchronization is functionally 
important [57]. 

Materials and Methods 



Synchronization stability analysis. The stability analysis can be carried out using 
a master stability approach based on linearizing Eq. [^around synchronous states 
1281 . We apply a coordinate transformation to reduce the Laplacian matrix L to 
its Jordan canonical form and decompose the variational equations into components 
along the corresponding eigenspaces fOl llOI . This leads to a master stability equation, 

y = [DY{s) - ai?H(s)]y, [9| 

wine re Cf, = sAj for the eigenspace corresponding to Ai . The stability function A((2) 
is the maximum Lyapunov exponent for the solution y{t) = of Eq. [9] The Aij 
are not required to be nonnegative and, since the derivation of A(tl) is based on 
a Jordan form, L is not required to be diagonalizable either [OlllOI . Geometrically, 
the stability condition is that al] the (possibiy comp]ex) numbers £A2, ■ • ■ , cAn 
lie in the region {o, ^ C : A(ci) <C 0} of the complex plane. Since A(ct) is 
generally positive on the left side of the imaginary axis, we consider only the networks 
whose Laplacian eigenvalues have nonnegative real part, which ensures that complete 
synchronization is possible. 

Dynamical robustness of optimal networks. For a given stability function A(a), 
a network satisfying condition[3]has the maximum rate of exponential convergence to 
synchronization among all networks. This is so under the assumption that there is a 
physical limit KI to the coupling strength of individual links, sAij ^ A/, and that 
the stability function A(a) is continuous and monotonically increasing with respect 
to the distance from the real axis, which appears to hold true for all known examples 
of A(a) in the literature [58] . From the limit on coupling strengths, it follows that 
the real part of s\i is bounded by a constant KIn{ri — 1), which, combined 
with the assumption on A(cf,), implies that the exponential rate of convergence is 
completely determined by the value of A(a) in the interval [0, Mn{n — 1)] 
on the real axis. In this interval A(a) has a global minimum at some a*, and 
thus r* := — A(ct*) > is the maximum possible rate of exponential conver- 
gence. If the network satisfies condition [3] all perturbation eigenmodes converge at 
the maximum rate when we choose £ = In contrast, if the network vio- 

lates condition!^ there must be at least one eigenvalue A^ such that sXi 7^ (2* 
[excluding the exceptional situation where multiple values of sXi fall precisely on 
multiple global minima of A(a)], resulting in an eigenmode that converges at a 
slower rate r := — A(£Ai) < r*. Therefore, although optimal networks may 
suffer from initially slower convergence to synchronization that is polynomial in time 
|?||10| . the long-term convergence rate is dominated by T and is faster than for any 
other network. Indeed, the deviation from the synchronous state can be written as 
P{t)e. ^ ^ for an optimal network and as Q{t)e, for a suboptimal network, 
where P{t) and Q{t) are polynomials. The ratio between the deviations in the two 
cases is then 

P{^^-{r*-r)t _^ Q as t 00 [10] 

and, in particular, is less than 1 for sufficiently large t, implying that the deviation 
will eventually become smaller for the optimal network. 

Laplacian spectrum of generalized complements. We show that if the eigen- 
values of L are 0, A2, ■ • • 5 An (counting multiplicity), then the eigenvalues of 
the Laplacian matrix of the generalized complement (defined by Eq. [TJ are 
0, no. — A2, . . . , no. — \n. This result follows from the relation 

fl(L^,x) = (—1)^"'"^ ~~ [11] 

na — X 

where x) = det(_L — x/) is the characteristic polynomial of the matrix 

L and / is the n X Tl identity matrix. We derive this relation by following the 
strategy of the proof of Lemma 2.3 in Ref. . 59 , for undirected networks with non- 
negative link strengths, to now consider directional and weighted links, possibly with 
negative strengths. From the definition of the complement transformation, we can 
write L + = TlCil — (xJ, where J is the n X 71 matrix with every entry 
equal to one. Using this and well-known properties of the determinant, we have 

//(L^, x) = det(L'^ — xl) 

= det(na/ — aJ — L — xl) 

= (-1)" det [L^ + a J - (na - x)l] 

= {-lY^(L^ + aJ.na - x), [12] 
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where L denotes the transpose of L. Eq. |11| then follows from the fact that 
^{L/^ + aj, z)/(na — z) = —fl(L^,z)/z whenever L has row sums 
equal to zero (which is the case, since L is the Laplacian matrix of a network). To 
prove this fact, we will use elementary row operations, which do not change the de- 
terminants. First, we replace the first row of the matrix L/^ + CtJ — zl by the 
sum of all rows, making each entry in this row TICX — Z. Next, we subtract this row 
multiplied by Ck/(^nct — 2) (which makes it a row with all entries equal to Ct) from 
the remaining rows, canceling out the contribution from the term CK</ in these rows. 
We denote the resulting matrix by AIi. Finally, we take the matrix L — zl and 
replace the first row by the sum of all rows. This results in a matrix ^2 with exactly 
the same entries as AIi except for the first row. For this matrix this is a row with all 
entries equal to — Z rather than TICX — Z. Dividing the first rows of Mi and AI2 
by TlOi — Z and — Z, respectively, which will scale the determinants accordingly, we 



see that 

^(L^ + aJ, z) _ det{L^ -\- aJ - zl) _ det(Afi) 
na — z no. — z no. — z 

_ det(M2) _ det(L^ - zl) _ fi{L^,z) 

—z —z z 

[13] 

As an immediate consequence of this result, the normalized standard deviation G 

changes to , under the complement transformation, while the total link 

° <y.nyi\ — \) — n\ ^ 

strength m is mapped to OLn(n — 1) — m. In particular, this implies that the 
complement of any optimal network is also optimal if Q >■ ^^^„2) ' 
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Supporting Information 



1. Laplacian spectrum of optimal networks 

Consider a network whose (possibly directional) links have in- 
teger strengths, and denote its Laplacian matrix by L. Here 
we show that if the network is optimal, i.e., the non-identically 
zero eigenvalues of L assume a common value A (Eq. 3 in the 
main text), then A must be an integer. This result will be 
valid even when the links are allowed to have negative integer 
strength. 

The characteristic polynomial of L can be written as 



det(I/ — xl) = —x{\ — x) 



+ --- + (-l)"a;". 



fSll 



where n is the number of nodes and I is the n x n identity 
matrix. Since L has integer entries, all the coefScients of the 
characteristic polynomial are integers, and hence A"~^ in the 
first term above is an integer. Denote this integer by k. Using 
m to denote the sum of all hnk strengths in the network, we 
have m = La = {n — 1)A, and hence k = [m/{n — 1)]"~^. 
Writing m/{n — 1) = s/t, where integers s and t do not have 



common factors, we obtain kt" 



Suppose p is a 



fact p is a factor of s. 
-1 



This implies that p" is a factor 



of kt^~^ = s"~^. Since s and t cannot have a common factor, 
p""^ must be a factor of k. Thus, any prime factor of k must 
actually appear with multiplicity n— 1, and hence we can write 
k = q"~^ where q is an integer. Therefore, A"~^ — k = q"~^ , 
and since A is real, we have A = g, an integer. 



2. Perturbation of Laplacian eigenvalues 

Suppose that the Laplacian matrix Lq of a given network of 
n nodes has an eigenvalue Ao 7^ with multiplicity k < n — 1. 
Consider a perturbation of the network structure in the form 
L — Lo + 5Li, where 5 is a small parameter and Li is any 
fixed Laplacian matrix representing the perturbed links. We 
do not need to assume that Lq and Li have nonnegative en- 
tries, making our result valid even in the presence of negative 
interactions. Denote the characteristic polynomial of L by 
f{x, S) := det(L — xl), where / is the n x n identity matrix. 
Since Ao is an eigenvalue of Lo with multiplicity k, we have 



f{x,0) = {x-Xorg{x), 



\S2] 



where 3 is a polynomial satisfying (/(Ao) 7^ 0. Denote by 
A = X{5) an eigenvalue of L that approaches Ao as 5 — > 0. 
Here we show that the change AA := A — Ao of the eigenvalue 
induced by the perturbation scales as 



AA . 



\S3] 



if the derivative of the characteristic polynomial with respect 
to the perturbation parameter evaluated at 5 = is nonzero: 



dl 
dS 



s^o 



\S4] 



Through the Jacobi's formula for the derivative of determi- 
nants, this condition can be expressed as 



tr[(io-AoJ)'=-'5'(io)ii] /O, 



fS5l 



where tr{A) denotes the trace of matrix A. We expect this 
condition to be satisfied for most networks and perturbations. 



For the optimal networks satisfying Eq. 3 in the main text, 
it can be shown that Eq. IS5I is violated if Lo is diagonal- 
izable, but Lo is actually known to be nondiagonalizable for 
the majority of these optimal networks [1]. The scaling [S3] 
shows that, for a fixed S, the more degenerate the eigenvalue 
(larger k), the larger the effect of the perturbation on that 
eigenvalue. In particular, if the original network is optimal 
with A having the maximum possible multiplicity n — 1, the 
effect of perturbation is the largest. This, however, is so be- 
cause the optimal networks have significantly smaller a than 
suboptimal networks (even those with just one more or one 
less link), and therefore the perturbations of the optimal net- 
works would still be more synchronizable in general than most 
suboptimal networks. 

From Eq. [S2]it follows that 



df 



Ox X — Aq Ox^ x^Xq 

s^o s=o 



dx'' 



s=o 



= 0, 



but 



9V 



x = Xo 
6=0 



= fc! • g{\o) / 0. 



fS6l 



[S7] 



Using this to expand f{x, 5) around x — Ao and (5 = up to 
the fcth order terms, and setting a; = A, we obtain 



fiKS) _ df 



+ 



1 d"! 



dS x=Xo k\ dx'' 
s=o 



X = Xq 

S=0 



+ 0(5), [S8] 



where 0{S) includes all higher order terms. From the charac- 
teristic equation /(A, S) = det(L — XI) = 0, the left hand side 
of Eq. ISSl is zero, so taking the limit 5 — >■ leads to 



lim 

S-fO 



(AA)* 



1 



df 



which implies the scaling [ 



g{Xo) dS 
Iwhen condition 



x^Xo 
S=0 



[S9] 

lis satisfied. 



3. Complexity of optimal networks 

Here we first describe a systematic method for increasing the 
number of nodes n in an optimal binary interaction network 
(Aij — 0,1), while keeping the network optimal. Given an 
optimal network with X = k (which must be an integer by 
the result in Section □ above) , we construct a new network by 
adding a new node and connecting any k existing nodes to the 
new node. As a result, the Laplacian matrix has the form 



L = 








Lo 









Ml • • • M„ 


k 



fsioi 



where Lo is the Laplacian matrix of the original network, each 
Ui is either or —1, and ui -|- • • • -I- m„ — —k. Since L is a 
block triangular matrix, its eigenvalue spectrum consists of 
the eigenvalues of Lo, which are 0, k, . . . , k, and an additional 
k, which comes from the last diagonal element. Thus, the new 
network is optimal with X — k. 

We can argue that the number of optimal binary inter- 
action networks grows combinatorially with n. To this end, 
we first consider C{n), the number of distinct Laplacian ma- 
trices corresponding to optimal networks with n nodes. For 
each optimal network with n nodes and X = k, the above con- 
struction gives (^) different Laplacian matrices correspond- 
ing to optimal networks with n + 1 nodes. Using the bound 
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(fc) ^ which is vahd for k = l,...,n — 1, we see that 
(7(71 + 1) > n • C{n), which imphes C{n) > (n — 1)! and 
gives a combinatorially growing lower bound for C(n). Since 
two different Laplacian matrices may represent isomorphically 
equivalent networks, C{n) is an overestimate of the number 
of optimal networks with n nodes. However, given the gross 
underestimate coming from (^) > n and the fact that we used 
only one out of potentially many possible schemes for adding 
a node while keeping the network optimal, we expect that the 
number of optimal binary interaction networks with n nodes 
also grows combinatorially with n. 



4. Optimality for networks of heterogeneous units 

Consider a network of coupled non-identical units whose dy- 
namics is governed by 

n 

[Sll] 

where t represents the discrete time and e = e/d is the global 
coupling strength normalized by the average coupling strength 
per node, d = ^ Yli Ylj^^i - The dynamics of unit i follows 
Xi(t -|- 1) = F{xi{t), Hi) in the absence of coupling with other 
units and is assumed to be one-dimensional for simplicity. 
Variation in the parameter fii represents the dynamical het- 
erogeneity of the network, which we measure by the standard 
deviation cr^ defined by cr^, :— ^ Sfif, where S/ii := fii — fi 
and p, := -^"^il-i-i- Here we choose the signal function to 
be H{x,h) — F{x,iJ.), which leads to a natural generaliza- 
tion of coupled map lattices [2] to arbitrary coupling topol- 
ogy. For example, for the one-dimensional periodic lattice in 
which each unit is coupled only to its two nearest neighbors 
with unit strength, system ISllI reduces to the well-studied 
system Xi{t + 1) = (1 - e)F{xi{t), fii) + ^F{xi^i{t), + 
^F{xi+i{t),fii+i). 

We consider a nearly synchronous state [3], in which 
the deviation of the states of individual units around their 
average is small, i.e., 5xi{i) := Xi — x{t) is small, where 
x(t) := ^'^iXi{t). Following the strategy used in Ref. j3j 
for continuous-time systems, we linearize Eq. ISlll around the 
average state x{t) and the average parameter fi to obtain the 
variational equation in the vector form 



5x(t + 1) = (/ - eL){at5^{t) -f bt<5/i). 



fS12l 



where &-x.(t) = {5xi{t), . . . , 5xn{t))'^ is the state deviation vec- 



tor, 5/1 = (5/11 



, (5/i„) is the parameter variation vector. 



and I is the n x n identity matrix. Matrix L is the modified 
Laplacian matrix defined by Lij = Lij — ^"^i^Lkj [3], and 
we denote [at,6t] :— [^{x{t), fl), ^{x(t), fl)j. As a result of 
the linearization, the deviation 5x(t) can possibly diverge as 
t (X even when the state space for the network dynamics 
is bounded. Notice that (1, . . . , 1)"^ is an eigenvector of the 
matrix I — EL associated with eigenvalue one. The component 
of the linearized dynamics parallel to this vector is irrelevant 
for synchronization stability, since 5x(t) by definition does 
not have this component. We thus remove this component, 
keeping all other components unchanged, by replacing / in 
Eq. IS12I with L* defined by L*j — Sij — 1/n, which leads to 




Fig. SI. Master synchronization error function f2(j4') defined on the space of 
matrices. The color indicates the values of Q,(^X) with X = L for the two- 
parameter family of networks given by L = CiLx.^ + C2{L(^^ ^ — 2) 
and the map F[x,fj.) = 2x -\- ^ mod 1. Shown in black are the level curves 
{^[X) = -^^^ - For fixed error tolerance i?tol heterogeneity cr^, the region 
defined by ^(X) < iEtoi represents the set of networks for which synchronization 
error is within -Etol- "^he red closed curve indicates the edge of synchronization 
stability defined by p{L* — X^ < . Each dashed line represents the path 

X = eh for a fixed h, which indicates the locations of synchronization transition 
(red dots) and, more generally, how synchronization error changes as e is varied (in- 
sets). The same representation for the path along the Ci-axis shows that ^{Xi) = 
when X = L*, corresponding to the point (ci,C2) = (1/3,0) for this example 
(indicated by a red plus symbol). 

Any component along (1, . . . , 1)"^ will immediately vanish un- 
der multiplication of the matrix L* — eL, whose properties 
govern the evolution of synchronization error. 

As a measure of synchronization error, we use the standard 
deviation iJx{t) defined by cj^{t) := i "Y^^ 5x1{t). For a fixed 
and e, we define the maximum asymptotic synchronization 
error to be 

^^(L) ~ max limsupcri:(t), [S14] 

{Mi} t— >oo 

where the maximum is taken over all possible combinations 
of Hi for the given a^. We can explicitly compute Q.^{L) by 
iterating Eq. IS13I which leads to 

Q,e{L) = OfiClieL), where 



Q.{X) :— limsup 

T-^oo 



llak+T-t bT-t{L' -X) 



fSlSl 



5x(t +1) = (L* - £L){at5^{t) + btSfi). 



[S131 



where ||-|| denotes the spectral norm for matrices. Notice that 
the synchronization error is a linear function of cr^, which is a 
consequence of using Eq. IS13I The function tl{X), whose ar- 
gument is a matrix X encoding the network structure, can be 
interpreted as a master synchronization error function since 
its functional form is determined entirely by the map F{x,fi) 
and the averaged trajectory x{t), and is independent of the 
network structure. A sufficient condition for fl{X) to be fi- 
nite is p{L* — X) < e^"^, where p(-) denotes the spectral ra- 
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max iVieXi) 

2<i<n 




Stable synchronization 

A(a) < 

Fig. S2. Synchronization error function for undirected networks. The condition for 
the stability of synchronization is that the numbers e\2 ; • ■ ■ i (the red dots on 
the a-axis) all lie in the interval on which A{a) < 0. When this condition is satisfied, 
the synchronization error is determined by the highest point among the corresponding 
points on the curve ri^(cl) (the red dots on the top curve). For illustration we used 
the curves for F{x,^) = 2x + fi mod 1, which are S7°(a) = | ^^^rrf | 
A(a) = ln2 + In |1 - a|. 



dius of matrices, or equivalently, the maximum of the abso- 
lute values of the eigenvalues. Here v is the Lyapunov ex- 
ponent of the averaged-parameter map F{x, p.) along the av- 
erage trajectory x{t), i.e., v ~ lim ^Yl't=o ^''^l^i^i^)' 

li X = eL, this condition reduces to the stability condition 
for complete synchronization of the corresponding identical 
units, namely, A{eXi) < for i — 2, . . . ,n, where the stability 
function in this case is A(a) = -|- In |1 — a|. For example, if 
F{x,fj,) — 2x + jj. mod 1, the sum in Eq. IS15I converges to 
Q{X) = ||(L* -X)(2X-2L* when p{L* - X) < 1/2. 

The set of networks with a given synchronization error 
tolerance Stoi defined by ile{L) < Etoi is represented by the 
region {X : Q.{X) < in the space of matrices. In Fig. ISII 

we illustrate this using the two-parameter family of networks 
defined by L = CiLks +C2(I/C3,i — ^03,2)1 where Lxa, ^03,1, 
and Lc3 2 ''■rs the Laplacian matrices of the fully connected 
network of 3 nodes and the two types of 3-cycles: 



2 


-1 






2 






-1 


2 




-1 








1 












^C3,i = 




[S161 



For this family of networks, we can show that for X = L, we 
have 



Cl{L) = lim sup 



[S171 



where A — 3ci -f i\/3c2. In particular, we have fl{L) = | 2J~lI 
if F{x,fj,) — 2x + jj, mod 1, and this is used as an illus- 
trative example in Fig. ISII As approaches zero and 
the dynamical units become less heterogeneous, the region 
{X : Cl{X) < -^^l increases in size and approaches the region 



of stable synchronization given by p{L* — X) < e"" — 1/2. 
For a given network, the change in the synchronization error 
0.£{L) with respect to e can be understood as the change in 
the value of the error function Q{X) as X — eL moves along 
a straight line determined by L. From Eq. IS15I we expect 
in general that Cl{X) is a monotonically increasing function 
of p{L* — X), and hence Qe{L) is expected to decrease to 
a minimum at some e = e* and increases monotonically for 
e > £*, or monotonically decrease in the entire range of e for 
which f2e(L) is finite. This is indeed the case for the example 
considered here, as illustrated by the insets in Fig. ISII Thus, 
in general we define Q{L) :— inf^ Qs{L) as a measure of syn- 
chronizability of a network, as it gives the lower limit on the 
asymptotic synchronization error for non-identical units. For 
undirected networks, for which L is symmetric and each A; is 
real, diagonalization of L with orthogonal eigenvectors can be 
used to show that 



Qe{L) = ai_iQ{eL) — max Q''{eXi), where 

2<i<n 



Cl" (a) = lim sup 

T->oo 



T /t-1 

I n "-k+T-t 1 bT-t(l - a)' 



[S181 



under the stability condition A{eXi) < for i = 2, . . . , n. For 
such networks, the synchronization error can be determined 
visually from the error function Q''{a), which is a function of 
real numbers. This is illustrated in Fig. IS2I using the exam- 
ple of F{x,ii) = 2x + fi mod 1, for which we can show that 
n"(a) = and A(a) = ln2 -Flu |1 - a|. 

We now show that the class of networks with zero syn- 
chronization error for arbitrary heterogeneity of the individ- 
ual units consists of those that are optimal [i.e., satisfies 
A2 = • ■ ■ = An = A > 0, Eq. 3 in the main text] and 
have diagonalizable Laplacian matrix. First, to show that 
any network with zero synchronization error satisfies these 
conditions, suppose that Q{L) = for a given network. That 
is, for some e, we have Sx{t) — > as t — ^ 00 for arbitrary 
/ii, . . . , /i„ with a given a^. Then, letting f — > 00 in Eq. [S13] 
we conclude that we either have lim bt — lim ^(x(t), n) = 

or (L* — eL)Sfi = for all possible Sfi. The former can hold 
only in exceptional cases, such as when x(t) converges to a 
fixed point at which ^ is zero. We thus assume a typical 
situation in which the latter holds. In this case, using the fact 
that the row sum of L is zero and that Spi — 0, we can 
show that L* — eL must be equal to the zero matrix, and hence 
L = jL* , which is diagonalizable with eigenvalues 0, 4, . . . , i. 
Since in general L and L have the same set of eigenvalues and 
L is diagonalizable iff L is diagonalizable, L is diagonalizable 
and satisfies A2 = • • • = An = A > 0. 

Conversely, suppose that the network satisfies A2 = • • ■ = 
An = A > and the Laplacian matrix is diagonalizable. It can 
be shown that eL = L* if we choose e — 1/A, and therefore 
5x(t) = according to Eq. IS13I but we can actually prove 
stronger statement without the linear approximation involved 
in Eq. IS13I From Theorem 6 in Ref. [1], each node j either 
has equal output link strength to all other nodes {Aij = bj ^ 
for all i 7^ j) or has no output at all {Aij = bj = for all i). 
This implies that the adjacency matrix satisfies Aij = bj for 
all i and j with i 7^ j, and we have bj = X. If we choose 
e = 1/A, then Eq. ISllI becomes 



^ n n 

x^{t + l) = ^Y.^jF{x,{t),p,) = ^^,F(x,(f),/i,), [S19] 
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Fig. S3. Change in the degree distributions when enhancing synchronization with 
negative directional interactions. Distributions before (blue) and after (red) assign- 
ing negative strengths in random scale-free networks: [A] in-degree distribution for 
"7 = 2.6, (B) out-degree distribution for 7 = 2.6, (C) in-degree distribution for 
7' = 5, and (D) out-degree distribution for '7 = 5. The in-degree distributions are 
plotted using logarithmic binning, and the absence of a symbol implies that we ob- 
served no node with its degree in the corresponding bin. The out-degree distributions 
are shown in linear scale, along with insets showing the positive part of the distri- 
butions in logarithmic scale. Note that there is a significant number of nodes with 
negative out-degree (green arrow). All plots are averaged over 20 network realizations 
with 1,000 nodes. 

with "^^j ™i ~ 1- Thus, the state of node i is determined by 
the weighted average of the signals from all the nodes that 
have output and, more importantly, it is independent of i for 
aU t > 1, implying that 5x(t) = for all t > 1. There- 
fore, the system synchronizes in one iteration with zero error, 
despite the presence of dynamical heterogeneity, and hence 
i}{L) = 0. This indicates that the largest Lyapunov exponent 
for the completely synchronous state is —00, which is analo- 
gous to super-stable fixed points and periodic orbits observed 
in maps. 

Since the best networks for synchronizing non-identical 
maps satisfy A2 = • • • = A„ = A > 0, they too must have 
a quantized number of links: m — k{n — 1). For every n 
and every k — 1, . . . ,n, there is exactly one binary network 
{Aij — 0, 1) that has m — k{n~l) links and is capable of com- 
plete synchronization for non-identical maps, including the di- 
rected star topology (fe = 1) and the fully connected network 
{k — n). Note also that the above argument does not require 
that bj > for all j (Theorem 6 in Ref. [1] remains valid with- 
out this requirement). This implies that complete synchro- 
nization is possible even for networks with negative interac- 
tions. In addition, the ability of a network to completely syn- 
chronize non-identical units, with or without negative interac- 
tions, is invariant under the generalized complement transfor- 
mation defined by Eq. 7 in the main text. To see this, suppose 
that for a given network we have A2 = • • • = An = A > and L 
is diagonalizable. By Theorem 6 in Ref. [J], we have Aij — bj 
for all i and j with i ^ j, where bj = X. Using the defini- 
tion of the complement transformation, we have Aij = ct — bj, 
and Y.j{a-bj) = na-\>Q-iia> Applying Theo- 

rem 6 in Ref. [1] again, we see that the complement satisfies 
the same property: A2 = • ■ ■ = A„ = na — A > and its Lapla- 




Fig. S4. improving synchronization with negative bidirectional interactions. {A 
and B) Degree distributions before and after we apply two methods of assigning 
negative interactions to random scale-free networks with 7 = 2.6 and 5: {A) 
degree-based method and {B) eigenvector-based method. All four plots are gener- 
ated using logarithmic binning and averaged over 20 network realizations with 1,000 
nodes. The absence of a symbol indicates that we observed no node with its degree 
in the corresponding bin. (C) Reduction in An / A2 as a function of 7' for the two 
methods. The error bars indicate the average and the standard deviation over 20 
network realizations with 1,000 nodes. The individual realizations are indicated by 
gray and blue dots. 

cian matrix is diagonalizable. Therefore, in addition to binary 
networks, there are many networks with negative interactions 
that are guaranteed to have zero synchronization error. 

5. Degree distribution before and after enhancing 
synchronization with negative directional interactions 

We describe the change in the in- and out-degree distributions 
of the network as negative strengths are assigned to directional 
links to enhance synchronization, following the algorithm pre- 
sented in the main text. The in- and out-degree of node i are 
defined as Y^^^^Aik and Y^^^^Aki, respectively. Figure [Sl 
shows the results for random scale-free networks with 7 = 2.6 
and 7 = 5. They clearly illustrate that the large in-degree 
of many nodes are compensated by the negative interactions, 
creating a sharp cut-off in the distribution (orange arrows in 
panels A and C). In contrast, the out-degree distributions 
remain essentially unchanged, having a power-law tail with 
the same exponent (insets in panels B and D). Note that the 
algorithm can create negative out-degree nodes, as indicated 
by the green arrow in panel B, but this has no significant ef- 
fect since the in-degree distribution is the main factor that 
determines the stability of synchronous states. 
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6. Enhancing synchronization with negative bidirec- 
tional interactions 

Here we show that assigning negative strength to bidirectional 
links can also enhance synchronization significantly. This is 
implemented using two different algorithms. 

The first method is fast and is based on node degrees, 
similarly to the algorithm used in the main text for assigning 
negative directional interactions. In order to create negative 
interactions preferentially between nodes of large degrees, we 
first order the bidirectional links according to the product of 
the degrees of the two nodes connected by each link, from high 
to low values. Going through all the links in this order, we 
change the strength of each bidirectional link from +1 to —1 if 
the degrees of the two adjacent nodes do not fall below a con- 
stant, chosen here arbitrarily to be 1.7 times the mean degree 
of the initial network. We applied this procedure to random 
scale-free networks with minimum degree 5, generated by the 
configuration model. Figure IS41 4 shows the degree distribu- 
tion before and after assigning negative interactions for the 
scaling exponent 7 — 2.6 and 5. In both cases, the degree 
distribution remains essentially scale-free with the same ex- 
ponents. Denoting by A2 and A„ the smallest and largest 
non-identically zero eigenvalues of the Laplacian matrix L, 
respectively, we measure synchronization enhancement by the 
relative decrease in the ratio \n/\2, a standard measure of 
synchronization widely adopted for undirected networks [4]. 
Figure [S4]C shows that our method does produce significant 
enhancement for 7 less than about 5 (bottom curve), and that 
the more heterogeneous the initial degree distribution (smaller 
7), the more effective the algorithm. 

The second algorithm is slower but more effective, and it 
is based on the observation that the largest eigenvalue A„ is 
typically the one responsible for poor synchronizability in the 
undirected networks with heterogeneous degree distribution 
considered here [S]- At each step of this algorithm, we use the 
first-order approximation [6l for each bidirectional link to 
estimate the change in A„ that would be caused by changing 
the strength of that link from +1 to —1. We then choose a 
link with the largest predicted reduction in A„ and make its 
strength —1. Repeating this until the fraction of links with 
negative strength reaches a prescribed threshold (chosen here 
arbitrarily to be 0.2), we obtain a sequence of candidate net- 
works for improved synchronizatiorQ. From these networks 
we choose one with the smallest ratio A„/A2 as the output of 
the algorithm. Figure [S4l B shows the change in the degree dis- 
tribution before and after applying this algorithm to random 
scale-free networks with minimum degree 5 for 7 = 2.6 and 5, 



generated by the configuration model. For both values of 7, we 
notice a significant drop in the fraction of high-degree nodes 
(green arrows) , accompanied by the appearance of nodes with 
degree less than the minimum degree of the initial network 
(blue arrow). This change appears to be responsible for the 
reduction of A„/A2 by as much as about 65%, as shown in 
Fig. IS4I C (top curve) . Note that the synchronization enhance- 
ment achieved by this method is consistently larger than the 
first method based on link degrees (bottom curve) . The effec- 
tiveness of this method depends on the fact that the Laplacian 
eigenvalues remain real when the network is kept symmetric, 
which would not generally hold true if negative strength were 
assigned to directional links. 

Supporting Video 

http : //www. youtube . com/watch?v=3dMIlYyxmbw 

Networks with best synchronization. The first half of the 
movie shows the structural changes in networks with best syn- 
chronization properties (smallest a possible for a given num- 
ber of links) as directional links are removed one by one. We 
always choose a link that keeps the synchronizability highest 
(i.e., keeps a smallest). The second half shows how a changes 
in the process, revealing in particular that link removal can 
counterintuitively enhance synchronization. The node layout 
at each step was computed using the Kodama-Kawai spring 
layout algorithm 8 . We used the implementation of the al- 
gorithm in the Boost C-I--I- library through the Matlab 
interface provided by MatlabBGL |10) . 
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